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Preface 


This  study  analyzed  the  relative  motion  between 
satellites  arranged  as  a  constellation  in  an  attempt  to 
quantify  the  magnitude  of  any  gravity  induced  relative 
motion.  It  was  felt  that  knowledge  of  this  type  of  motion 
could  be  used  by  designers  of  satellite  constellations 
for  sizing  of  propulsion  systems  for  station-keeping 
and/or  for  deciding  what  orbital  parameters  to  use  for  the 
satellites  within  the  constellation. 

Motion  was  analyzed  between  two  satellites  in 
circular  orbits  at  altitudes  ranging  from  300  km  to 
6000  km  and  orbit  inclinations  from  40°  to  80°.  The  J2 
gravitational  harmonic  was  used  as  the  perturbing  force, 
and  results  from  18  different  test  cases  showed  that  rela¬ 
tive  motion  does  occur  as  an  undamped  oscillation  in  range. 
However,  the  magnitude  of  this  oscillation  under  J2  per¬ 
turbations  is  small  enough  that  the  relative  motion  need 
not  be  applied  as  a  constraint  on  satellite  constellation 
design.  The  method  employed  does  provide  a  means  of 
analyzing  other  perturbing  forces,  and  this  could  be  done 
to  further  substantiate  the  findings  of  this  study. 

I  wish  to  extend  my  thanks  to  Captain  Keith  Jenkins 
from  the  Systems  Concepts  Group  at  AFWAL  for  identifying 
this  topic  as  an  area  which  needed  investigation.  I  also 
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want  to  thank  Lt  Col  Joe  Widhalm  for  agreeing  to  be  the 
advisor  for  this  thesis  even  though  I  wasn't  a  "real 
astro  guy." 

—  Mark  J.  Buechter 
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Abstract 

Current  plans  for  systems  to  be  used  for  ballistic 
missile  defense  sometimes  call  for  using  satellites  that 
are  placed  so  as  to  form  a  constellation  that  can  continu¬ 
ously  monitor  specified  areas  of  the  Earth's  surface  or 
direct  weapons  against  attacking  missiles  or  warheads. 

This  study  analyzes  the  relative  motion  between  satellites 
within  such  a  constellation  under  gravitational  perturba¬ 
tions  caused  by  the  Earth's  equatorial  bulge  (oblateness). 
Relative  motion  is  calculated  using  a  system  of  equations 
which  describes  the  variation  of  relative  orbital  elements 
between  two  satellites.  The  cases  studied  simulate  the 
position  of  two  satellites  that  are  located  within  a  con¬ 
stellation  containing  ten  orbi'ts  with  ten  satellites  in 
each  orbit.  The  orbits  investigated  were  all  circular  with 
altitudes  ranging  from  300  km  to  1000  km  and  inclinations 
ranging  from  40  degrees  to  80  degrees.  Range  between  the 
satellites  was  oscillatory  with  deviations  from  the  average 
range  of  up  to  5  km.  The  results  vary  only  slightly  with 
changes  in  orbital  inclination  or  altitude.  These  results 
show  that  the  relative  motion  between  satellites  in  a  low 
altitude  constellation  caused  by  the  Earth's  oblateness 
does  not  significantly  affect  the  initial  geometry  of  the 


constellation. 


AN  ANALYSIS  OF  GEOPOTENTIAL  EFFECTS 


ON  SATELLITE  CONSTELLATIONS 


Background 


President  Reagan,  in  a  speech  given  in  March  of 
1983,  proposed  that  the  scientific  community  of  the  United 
States  help  develop  a  capability  for  ballistic  missile 
defense  (BMD)  to  render  nuclear  weapons  "impotent  and 
obsolete"  (3:39).  This  proposal  touched  off  debate 
throughout  the  country's  scientific  community  as  to  how 
technologically  achievable  or  economically  feasible  this 
type  of  defense  system  would  be. 

Both  sides  in  this  argument  have  put  forth  pro¬ 
posals  on  how  such  a  system  for  BMD  should  be  structured. 
All  of  these  options  call  for  some  type  of  airborne  or 
spaceborne  sensor  system  to  detect  incoming  missiles/war¬ 
heads,  but  the  various  options  differ  on  how  much  of  the 
actual  weaponry  should  be  placed  in  orbit. 

When  orbiting  weapons  are  needed  they  must  be 
placed  in  an  arrangement  that  allows  continuous  coverage 
of  those  areas  of  the  Earth's  surface  over  which  a  ballis¬ 
tic  missile  attack  might  occur.  When  several  satellites 
are  placed  into  orbits  so  that  they  can  provide  comple¬ 
mentary  capabilities  for  sensor  coverage,  data  relay  or 


any  other  purpose,  the  satellites'  position  relative  to 
one  another  will  need  to  be  specified.  Such  arrangements 
of  satellites  are  commonly  referred  to  as  constellations 
and  generally  consist  of  several  orbits,  all  of  equal 
geometry  and  inclination,  with  their  lines  of  nodes  spaced 
equally  around  the  Earth's  equatorial  plane.  Each  orbit 
contains  several  satellites  evenly  spaced  throughout. 
Figure  1  illustrates  such  a  constellation  of  satellites. 


Fig.  1.  Illustration  of  a  Satellite  Constellation 

The  satellites  may  be  the  actual  weapons  or  in  some  cases 
they  are  mirrors  used  to  redirect  a  laser  beam  emanating 
from  a  surface  based  weapon.  In  either  case,  the  relative 
positions  of  the  satellites  within  the  constellation  will 
determine  the  coverage  the  weapon  has  over  the  Earth’s 
surface.  Should  the  satellites'  spacing  within  the 


constellation  become  disrupted  there  is  the  possibility 
that  a  gap  may  occur  in  the  weapons'  coverage.  This  gap 
would,  of  course,  be  undesirable. 

It  appears,  however,  that  none  of  the  designers 
of  these  constellations  have  looked  at  how  perturbing 
forces  will  affect  the  constellations.  Since  all  of  the 
current  proposals  put  the  satellites  at  relatively  low 
altitudes,  the  orbits  are  considered  near-earth  and  will 
experience  perturbations  due  to  the  Earth's  non- spher ica 1 
shape.  These  perturbations  could  result  in  oscillations 
of  the  relative  positions  of  the  satellites  within  the  con¬ 
stellation  that  may  be  inconsequential  or  perhaps  will 
lead  to  total  disruption  of  the  constellation. 

Research  Problem 

It  is  necessary  to  know  how  the  actual  gravita¬ 
tional  field  of  the  Earth  will  affect  the  movement  of  satel¬ 
lites  within  a  constellation  to  determine  if  there  are 
adverse  changes  in  the  weapons  coverage  of  the  total  system. 
There  may  be  some  altitudes,  inclinations  or  combinations 
that  need  to  be  avoided  because  they  tend  to  accentuate 
the  relative  motion  between  the  satellites,  or  that  can 
be  capitalized  upon  to  lessen  any  relative  motion  problem. 

If  some  constellation  geometries  are  deemed  vital  for  the 
weapon's  coverage  but  run  contrary  to  the  satellites' 
natural  motion,  appropriate  propulsion  systems  will  need 
to  be  included  in  the  satellites'  design. 
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Previous  Research 

A  literature  search  for  previous  research  con¬ 
cerning  satellite  motion  within  a  satellite  constellation 
revealed  nothing  directly  addressing  this  problem.  How¬ 
ever,  this  initial  screening  did  reveal  work  on  the  motion 
of  one  satellite  that  is  very  close  to  another  (7;  8) . 

This  work  was  directed  at  applications  such  as  the  proposed 
free-flying  microgravity  lab  to  be  used  with  the  U.S.  space 
station,  but  also  appeared  applicable  to  the  satellite  con¬ 
stellation  problem.  The  foundation  of  the  analysis  was 
rendezvous  theory,  and  since  the  rendezvous  problem  also 
involves  solving  for  the  relative  motion  between  two  orbit¬ 
ing  objects,  research  into  the  solution  methods  for  the 
rendezvous  problem  seemed  an  appropriate  starting  point 
for  the  constellation  problem.  Early  studies  of  the 
rendezvous  problem  give  approximate  solutions  and  deal 
only  with  special  conditions,  but  the  more  recent  work 
gives  improvements  on  these  techniques. 

One  of  the  earliest  attempts  at  solving  the  rela¬ 
tive  motion  problem  inherent  in  the  rendezvous  of  two 
orbiting  satellites  was  accomplished  by  Eggleston.  In  his 
study,  Eggleston  looked  at  how  to  calculate  the  minimum 
change  in  velocity  needed  to  achieve  rendezvous  between 
two  satellites,  and  listed  several  sets  of  equations  which 
can  be  used  to  solve  this  problem  depending  upon  the  amount 


of  accuracy  needed  (4:27).  However,  this  investigation 
was  limited  to  circular  orbits  about  a  spherical  earth 
(4:10)  . 

Anthony  and  Sasaki  further  developed  the  solution 
for  the  rendezvous  problem  by  allowing  the  use  of  orbits 
with  small  eccentricities  and  by  including  second-order 
terms  in  the  equations  (1:1669).  These  additions  to  the 
solution  method  allowed  the  problem  to  be  solved  for  either 
large  or  small  relative  distances  (1:1673). 

The  research  by  Anthony  and  Sasaki,  and  Eggleston 
provided  only  approximate  solutions  to  the  relative  posi¬ 
tion  problem,  but  Lancaster  (5:1878)  showed  how  to  solve 
the  problem  exactly.  Like  the  previously  discussed  works, 
however,  Lancaster's  solution  was  also  based  on  a  spherical 
earth  assumption. 

All  of  these  methods  for  solving  the  relative  motion 
problem  were  based  on  assumptions  that  simplified  the  solu¬ 
tion  method.  The  one  common  assumption  was  that  the  satel¬ 
lite  was  orbiting  a  spherical  planet,  and  therefore  the 
gravitational  field  of  that  planet  was  uniform.  Since  the 
Earth  has  many  nonuniformities  in  its  mass  distribution, 
the  primary  one  being  a  bulging  at  the  equator,  an  object 
in  orbit  experiences  perturbing  forces  as  it  passes  over 
different  areas  of  the  Earth's  surface.  To  solve  for  the 
true  motion  of  the  satellite  when  these  forces  are  added 
to  the  satellite's  equations  of  motion  requires  the  use  of 
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a  computer  for  numerical  integration.  Unfortunately,  when 
the  two  satellites  of  interest  come  close  together,  and 
their  position  vectors  become  almost  equal,  much  precision 
in  the  solution  is  lost  because  of  the  requirement  to 
subtract  these  nearly  equal  vectors. 

Nacozy  and  Szebehely  developed  a  technique  for 
finding  relative  position  and  velocity  using  numerical 
integration  that  retains  the  precision  of  the  solution  when 
two  nearby  satellites  are  studied  (7:449).  They  used  the 
basic  technique  outlined  by  Eggleston,  but  improved  it 
to  allow  for  non-circular  orbits  (7:449).  Nacozy  and 
Szebehely' s  study  was  done  in  conjunction  with  a  study  by 
McKenzie  in  which  eight  different  sets  of  equations  were 
evaluated  for  computational  efficiency  (6:10).  Efficiency 
was  measured  by  the  number  of  computations  done  by  the 
computer  during  numerical  solution  of  the  equations.  The 
equations  identified  in  McKenzie's  report  as  Set  6  apply 
the  Nacozy-Szebehely  technique  to  the  classic  method  of 
Cowell  for  perturbations  and  are  applicable  to  the  case 
where  the  two  satellites  are  in  nearly  equal  orbits  (6:15). 

Van  der  Ha's  work,  the  foundation  of  which  can  be 
traced  to  all  of  the  previously  cited  works,  developed 
equations  for  the  motion  between  two  satellites  requiring 
no  simplifying  assumptions.  Van  der  Ha's  technique  trans¬ 
forms  the  satellites'  position  and  velocity  vectors  into 
six  elements  which  completely  describe  the  relative 
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1  Scope  of  the  Research 

h  Because  the  current 

proposals  for  satellite  con- 

-  stellations  to  be  used  for  b 

allistic  missile  defense  (BMD) 

M 

a  call  for  low  altitude  orbits 

(less  than  one  earth  radius) , 

irregularities  in  the  Earth's  gravitational  field  will  be 
a  major  perturbing  force.  Therefore,  this  investigation 
focused  on  gravitational  perturbations  and,  more  specifi¬ 
cally,  was  limited  to  inclusion  of  the  primary  oblateness 
term  (the  ^  harmonic)  from  the  mathematical  expansion  of 
the  geopotential. 

Captain  Keith  Jenkins  from  the  Air  Force  Wright 
Aeronautical  Laboratories,  who  originally  posed  this 
research  question,  noted  that  a  constellation  of  10  circular 
orbits  with  10  satellites  per  orbit  at  an  altitude  of 
1000  km  is  being  used  as  a  baseline  by  many  BMD  designers. 
This  constellation  geometry  was  used  as  a  baseline  from 
which  several  test  cases  were  drawn.  The  basic  investiga¬ 
tion  looked  at  the  relative  motion  between  two  adjacent 
satellites  in  the  same  orbit,  and  at  the  motion  between 
two  satellites  in  different,  but  adjacent,  orbits.  Then 
the  altitude  and  inclination  of  the  orbits  were  varied  to 
determine  the  sensitivity  of  the  results  to  these  param- 


II.  Methodolo 


Discussion  of  Method  Used 

The  method  chosen  to  conduct  this  study  of  relative 
motion  between  satellites  in  a  constellation  was  to 
numerically  integrate  a  set  of  equations  of  motion  for  one 
satellite  and  a  set  of  equations  describing  the  change  in 
the  state  of  the  second  satellite  with  respect  to  the  first. 
This  method  is  described  by  Van  der  Ha  (8)  and  major  points 
will  be  described  here.  The  core  of  this  method  is  the 
development  of  a  set  of  difference  elements  between  the  two 
satellites,  and  equations  for  the  time  rate  of  change  of 
these  difference  elements.  Difference  elements  are  used 
to  avoid  losing  precision  in  the  solution  caused  by  the 
subtraction  of  two  nearly  equal  numbers  (8:293).  This 
methodology  allows  analysis  of  relative  motion  regardless 
of  how  close  the  two  satellites  approach  one  another. 
Therefore,  this  methodology  could  be  applied  to  satellite 
constellations  containing  large  numbers  of  satellites  as 
has  been  suggested  by  the  distributed  assets,  or  "swarm," 
concept  for  deployment  of  space-based  ballistic  missile 
defenses . 

To  allow  for  circular  and/or  equatorial  orbits, 
the  state  of  the  satellites  is  described  using  a  set  of 


orbital  elements  different  from  the  classical  elements. 


These  elements  are  listed  as  Equation  (1) . 

"a  =  (h ,  f , g,  j  ,k,L)  (1) 

These  elements  are  written  in  terms  of  the  classical 
elements  as 

h  =  /uT 

f  =  e  cos  6 
g  =  e  sin  c 
j  =  tan  (i/2)  cos  f. 
k  =  tan  (i/2  )s inf- 

L  =  6  +  u  +  C  (2) 


where  i  is  the  semi-lattus  rectum,  e  is  the  eccentricity, 
i  is  the  orbit's  inclination,  '  is  the  longitude  of  the 
ascending  node,  ...  is  the  argument  of  perigee,  and  -  is  the 
true  anomaly  of  the  satellite. 

Conversely,  the  classical  elements  can  be  written 
in  terms  of  the  elements  in  Equation  (1)  as 


e  =  Ji2  +  g2 
-  =  arctan(g/f) 

1  =  2  arctan(  J j2  +  k2) 

2  =  arctan(k/j) 

w  =  l  -  e  -  n  (3) 


10 


s. 


vV- 


Aa“  =  (Ah,Af  ,Ag,A  j  ,  Ak  ,  AL)  (4) 

These  difference  elements  are  defined  as  the  dif¬ 
ference  between  the  orbital  elements  of  the  two  satellites 
where  one  of  the  satellites  is  designated  as  the  reference. 
The  relationship  between  the  positions  of  the  two  satel¬ 
lites  is  shown  in  Figure  2  where  the  coordinate  system 

is  an  inertial  frame  centered  at  the  Earth's  mass 
center  and  the  X^  and  X 2  axes  are  in  the  equatorial  plane. 


Fig.  2.  Relationship  Between  Position 
Vectors  of  Two  Satellites 
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Therefore  one  can  write 


+  A  r 


(5) 


or 


(6) 


This  relationship  can  also  be  applied  to  the  ele¬ 
ment  sets  so  that 


and  therefore 

Ah  =  h2  -  h1 


(7) 


(8) 


A  set  of  six  differential  equations  can  be  written 
describing  the  change  in  the  orbital  elements  with  time 
which  include  the  components  of  the  perturbing  force  on 
the  satellite.  These  force  components  are  with  respect 
to  a  reference  frame  centered  at  the  satellite  with  the 
x-axis  along  the  position  vector  from  the  Earth  to  the 
satellite,  the  z-axis  along  the  orbit  normal,  and  the 
y-axis  in  the  direction  which  completes  the  right-handed 
triad  as  shown  in  Figure  3. 


Fig.  3.  Satellite  Centered  Reference  Frame 


Development  of  these  equations  is  through  the 
technique  of  variation  of  parameters,  and  the  six  equations 
are  as  follows  (8:292) 

h  =  rFy 

f  = -hg/r2  +  2hFy/y 

g  =  hf/r2  +  hF  /y  +  grF  /h 
x  y 

j  =  rF  cos L  (l+j2+k2) /  (2h) 
z 

k  =  rF  sin  L  (l+j2+k^) / (2h) 
z  J 

L  =  h/r2  +  rFz  ( j  sin  L-k  cos  L)  /h  (9) 

where 

r  =  h2/  [y  (1+f)  ] 


and  F  ,  F  and  F  are  the  components  of  the  perturbing 
x  y  z 

force  along  the  corresponding  local  axes  shown  in  Figure  3 
(8:292)  . 


A  partial  application  of  these  rules  to  Equa¬ 
tion  (10)  leads  to 

Af  »  -(  (g+Ag)A  (h/r2)  +  hAg/r2]  +  2A (hF  ) /u  (12) 

where  the  first  term  in  Equation  (10)  has  been  grouped  so 
2 

that  x  =  h/r  and  y  =  g  for  use  with  Equations  (11). 

Similar  applications  of  these  rules  allow  one  to 
derive  the  following  full-precision  equations  for  the 
variation  of  all  of  the  difference  elements  (C:293) 


Ah  «  A(rFy) 

Af  =  - [ (g+Ag) A  (h/r2)  +  hAg/r2]  +  2A  (hF  ) /u 
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Ag  =  [  (f+Af  )A  (h/r  )  +  hA  f /r  ]  +  A  (hF  ) /u 

X 

+  [  ( g+A  g )  A  Gy  +  GyA  g  ] 


Aj  =  h  [  (1+J+AJ)AG  +  G  AJ]  cos  (L+AL) 

+  H  G  { 1+J )  A  ( cor  L) 
z 


Ak  =  k  [  ( 1+J+A J) A G  +  G  AJ]  sin  (L+AL) 


+  h  Gz  (L+J  )  A  (sin  L) 


AL  =  A  (h/r  )  +  (jsinL-kcosL)AG  +  (G+AG)[Ajsin  (L+l  L) 

z  z  z 

-  Ak  cos  (L+AL)  +  jA(sinL)  -  kA  (cos  L)  ] 


(13) 


where 


2  2 
J  =  J  +  K 


Gj  =  rF^/h,  j=y , z 
AJ  =  A  ( j2)  +  A  (k2) 


AG.  =  [A  (rF .  )  -  G.Ah]  /  (h+Ah)  ,  i=y,z 

X  XX 

Ar  =  [A  (h2 )  -  h2A  f  /  ( 1  +  f )  ]  /  [„  (1  +  f+A  f ) 


(h/r2 )  =  [Ah-hA  (r2) /r2]  /  (r+Ar)  2 


It  should  be  noted  that  the  difference  force  com¬ 
ponents  AF^,  AFy  and  AF^  should  not  be  interpreted  as  the 
vector  difference  of  the  two  perturbing  forces,  -  F^, 
written  in  the  local  reference  frame  of  Satellite  1. 
Rather,  each  difference  force  component  is  simply  the 
scalar  difference  of  the  corresponding  components  of  each 


vector  written  in  its  own  local  frame  (8:292).  For 
example,  if  the  local  unit  vectors  for  Satellite  1  are 
x,y,z  and  the  local  unit  vectors  for  Satellite  2  are  i,j,k, 
the  local  perturbing  forces  would  be  written  as 

F1=Xx+Yy+Z  z 

Eh,  =  Ii  +  Jj+Kk  (14) 

The  difference  force  components  are  then 

F  =  I  -  X 
x 

F  =  J  -  Y 

y 

Fz  =  K  -  Z  (15) 

The  perturbing  forces  themselves  are  derived  from 
an  equation  for  the  potential  energy  of  the  satellite. 

This  potential  function  follows  from  Poisson's  equation 
and  can  be  written  in  spherical  coordinates  as 


V(r,  i  ,  :)  = 


™  n 

— —  -  ( R  /r ) n  P ™  (cos 

n=0  m=0  e 


x  [S  sin (m  i)  +  C  cos  (mt)  ] 
nm  nm 


(16) 


where  r  is  the  length  of  the  position  vector  from  the 
center  of  mass  of  the  Earth  to  the  satellite,  i  is  the 


co-latitude  and  i  is  the  longitude  as  shown  in  Figure  4 


Fig.  4.  Definition  of  Co-latitude  and  Longitude 

Other  constants  appearing  in  the  geopotential 

function  are  the  universal  gravitational  constant,  G,  the 

mass  of  the  Earth,  M,  and  the  Earth's  mean  equitorial 

radius,  R  . 

e 

The  perturbing  force  is  then  calculated  as 

F  =  - VV 

The  potential  function  can  be  written  in  terms 
of  the  inertial  cartesian  coordinates  by  observing  from 
Figure  4  that  cosY  =  X^/r.  Thus  Equation  (16)  can  be 
written 

V  =  -3MRe2  J2[l/3  -  (X3/r)2)  /  (2r3)  (17 


where  M 


GM. 


■  v»  y  jr*jr»TT»  y^r*yT-yjr»ini  ir»jr-jT»jr»ir»jr»r^.y^j 


Solving  for  force  yields 

F*  =  -TV  =  -£  [1-5  (X3/r)  2]  r/r5  -  2-:X3U3/r5  (18) 

2 

where  £  =  3nj  R  / 2. 

2  e 

However,  the  quantities  F  ,  F  and  F  needed  in 

x  y  z 

Equations  (5)  are  components  of  the  force  vector  written 
in  the  local  frame,  so  the  force  vector  must  be  transformed 
into  the  local  reference  frame.  This  transformation  can 
be  written  as 


U, 

X 

1 

u 

y 

=  [A] 

U2 

u 

U-, 

z 

3 

_  _ 

_  _ 

where  th  are  the  components  of  a  vector  in  the  inertial 
frame  and  u^  are  the  local  frame  components.  The  com¬ 
ponents  of  the  transformation  matrix  are  (8:306) 


"ll 

— 

kn^3  +  cos  L 

1 12 

= 

“j"l3  +  s;'-nL 

"l3 

= 

2  ( j  sin  L  -  k  cos  L)  /  (1+J) 

"21 

= 

ka^3  -  sin  L 

"22 

= 

-jot  +  cos  L 

a23 

= 

2  ( j  cos  L+k  sin  L)  /  (1+J) 

"31 

= 

2k (1+J) 
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32 


=  -2 j / (1+J) 


i33  =  (1-J) / (1+J) 


(20 


Performing  this  transformation  yields  (8:298) 


Fx  -  -i<l-3:.132,/r4 


Fy  ~  2;a13i23/r 


Fz  =  -2^i3:‘33/r 


(21 


The  difference  components  are  found  using  the 
relationship 


JFi  =  <F2)  -  (F^  ,  i  =  x  ,y ,  z  (22 

i  i 

and  the  relationships  of  Equations  (11). 

!Fx  =  --[3i(132)  +  (  1-3  ,1  32)  1  (r4)  /r4]  /  (r+ir)  4 

JFy  =  2  [  li3->23“  (r  *  /r  “  (  ‘13  <23)  ^  /  <r+-r> 

iFz  =  2  •;  [  .  13  ,33i  (r4  )  /r4  -  1  (  ;  ]_  3  '  3  3 )  1  /  ( r  +  1.  r )  4  (23 


Imp lemen ta tion  of  Equation s 

To  solve  the  equations  for  the  variation  of  the 
orbital  elements  and  the  difference  elements  a  computer 
program  was  written  using  Fortran  77  and  a  numerical  inte 
gration  package.  The  numerical  integrator  employed  was  a 
fourth  order  predictor-corrector  called  HAMING  which  was 


provided  by  Dr.  William  Wiesel  of  the  Air  Force  Institute 
of  Technology,  Wr ight-Patter son  Air  Force  Base,  Ohio. 

This  program  provided  a  time  history  of  the  orbital  ele¬ 
ments  of  Satellite  1  and  of  the  difference  elements,  and 
was  verified  using  two  test  cases  whose  results  could  be 
predetermined.  Both  verification  tests  involved  unper¬ 
turbed  motion. 

The  first  test  case  used  two  satellites  in  the  same 
orbit  but  separated  by  36°  of  true  anomaly.  The  orbital 
parameters  in  terms  of  classifical  elements  were 

1  =  6778  km  ,  =  0° 

e  =  0  n  =  0° 

i  =  45°  ei  =  20° 


where  subscripts  denote  parameters  for  the  first  and 
second  satellites.  For  unperturbed  motion,  the  separation 
between  the  satellites  should  remain  constant  because  the 
satellites  are  in  circular  orbits.  This  was  noted  by 
observing  the  distance  between  the  satellites  at  each  time 
step  during  the  integration.  The  initial  separation  dis¬ 
tance  was  calculated  using  the  equation 

R  =  2  J  sin  [(c2-ei)/2]  =  4189.034  km  (24) 

where  R  is  the  range  between  the  two  satellites. 


When  this  first  test  case  was  run  on  the  computer 
program,  the  distance  between  the  satellites  did  remain 
constant  at  the  predicted  value  of  4189.034  km. 

The  second  verification  test  case  consisted  of 
having  two  satellites  in  circular  orbits  of  equal  altitude 
but  at  different  inclinations.  The  two  orbits  intersected 
along  their  lines  of  nodes  and  the  test  run  was  begun  with 
both  satellites  together  at  the  ascending  node.  The 
initial  conditions  were 

1  =  6778  km  *  =  0° 

•  e  =  0  ft  =  0 0 

11  =30°  6  =  0° 

12  =  60° 

For  unperturbed  motion  the  range  between  the  two 
satellites  will  oscillate  from  the  initial  condition  of 
zero  to  a  maximum  value  after  one-fourth  of  a  period  and 
back  to  zero  after  one-half  period.  This  sequence  will 
repeat  every  half  period.  The  relationship  between  the 
satellites  after  one  quarter  period  is  illustrated  in 
Figure  5. 

The  range  at  this  time  is 

R  =  2  Zsin[ (i2-i1) /2]  =  3508.551  km  (25) 

The  time  history  of  this  test  case  is  shown  in  Figure  6. 


Since  there  was  no  real  data  against  which  the 

computer  program's  output  could  be  compared  for  perturbed 

relative  motion,  the  author  attempted  to  verify  Van  der  Ha's 

force  equations  by  deriving  the  perturbing  force  components 

expressed  in  the  satellite's  local  frame.  However,  the 

force  components  as  expressed  by  Equations  (21)  could  not 

be  obtained,  so  intermediate  solutions  for  F  ,  F  ,  and  F 

x  y  z 

were  implemented  in  the  computer  code. 

Several  cases  were  run  using  this  new  coding  for 
the  perturbed  force  components  and  the  original  version  of 
the  computer  program  which  used  the  expressions  shown  in 
Equations  (21) .  Identical  results  were  obtained  from 
both  versions  of  the  program,  so  Equations  (21)  were  con¬ 
sidered  to  be  accurate  and  coded  properly. 

Since  the  results  of  the  verification  runs  were 
as  anticipated  and  the  code  appeared  to  contain  no  syntax 
errors,  the  model  was  considered  to  be  verified  and  ready 


III.  Discussion  of  Results 


With  the  equations  programmed  and  the  computer 
routine  verified  for  proper  execution,  a  series  of  test 
cases  was  developed  to  examine  the  relative  motion  problem 
and  its  sensitivity  to  changes  in  orbit  altitude  and 
inclination.  As  mentioned  in  the  opening  chapter,  two 
basic  scenarios  were  used  to  evaluate  relative  motion 
within  a  constellation:  two  satellites  in  the  same  orbit 
but  separated  by  some  amount  of  true  anomally,  and  two 
satellites  at  identical  positions  within  two  different 
orbits.  In  the  second  scenario  the  orbits  were  of  identi¬ 
cal  geometry  but  had  different  ascending  nodes.  While  the 
second  scenario  geometry  could  result  in  satellites  being 
in  orbits  that  would  result  in  a  collision  (polar  orbits) , 
this  situation  was  not  investigated.  The  baseline  orbit 
investigated  was  circular  with  a  radius  of  7378  km  and 
inclined  60  degrees  to  the  equatorial  plane. 

First  Scenario 

For  the  first  scenario  the  two  satellites  were 
placed  36  degrees  apart  within  the  orbit  and  the  relative 
motion  calculated  for  ten  orbits.  The  difference  in  true 
anomally  of  36  degrees  between  the  satellites  corresponds 
to  the  spacing  between  te..  satellites  equally  spaced  within 
one  orbit.  The  results  of  this  test  case  using  the 
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baseline  orbit  are  shown  in  Figure  7,  where  range  is  simply 
the  magnitude  of  the  relative  position  vector  between  the 
two  satellites. 

Using  Equation  (24)  the  unperturbed  range  was 
calculated  to  be 

R  =  2  *  sin(  (Sj-Oj) /2)  =  4560  km  (26) 

However,  as  shown  in  Figure  7  the  average  value  around 
which  the  range  oscillates  remains  constant  at  the  value 
of  4624.5  km.  The  difference  between  these  two  values  can 
be  explained  in  part  by  the  way  in  which  the  initial 
conditions  for  the  computer  program  were  calculated  and 
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partly  by  the  fact  that  the  perturbed  orbit  will  be 
slightly  non-circular. 

The  initial  conditions  used  as  input  to  the  main 

computer  program  were  established  by  first  defining  the 

initial  orbit  and  the  satellite's  position  within  it  using 

the  classical  orbital  elements  (;.  ,e  ,  i  ,  ,  ,  -  )  .  These  orbital 

elements,  transformed  into  the  element  set  and  defined  in 

Equation  (1)  and  the  equations  for  the  variation  of  the 

elements.  Equations  (13) ,  were  then  integrated  for  one- 

tenth  of  a  period.  The  values  for  the  orbital  elements  at 

this  time  were  used  as  the  initial  conditions  for  the  first 

satellite.  Integration  for  another  one-tenth  of  a  period 

yielded  another  set  of  orbital  elements  to  use  as  initial 

conditions  for  the  second  satellite.  This  was  done  to 

ensure  that  the  initial  orbital  elements  for  both  satellites 

were  in  the  same  perturbed  orbit.  A  vector  of  these  orbital 

elements,  (an)  and  (a~J  <  was  then  used  as  the  initial 
1  U  2 

state  vector  for  the  main  computer  program. 

However,  it  was  not  possible  to  obtain  al  and  "a" 

U1  U2 

exactly  one-tenth  of  a  period  apart  because  of  the  inte¬ 
gration  step  size.  For  example,  the  total  period  for  an 
unperturbed  circular  satellite  is  found  using  the  equation 

TP  =  (2-/  ,T )  r1’5  (27) 

where  r  is  the  radius  of  a  circular  orbit.  For  an  orbit 
with  a  radius  of  7378  km  the  period  will  be  6307  seconds. 
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If  one-tenth  of  this  time  is  rounded  off  to  630  seconds, 
there  is  a  round-off  error  of  .7  seconds  which  can  be 
transformed  into  an  angular  error  by  multiplying  by  the 
angular  velocity  of  the  satellite.  Angular  velocity  is 
calculated  as  follows: 

^SAT  =  360°/6307  sec  =  -05708  deg/sec  (23) 

Multiplying  ^  by  the  round-off  time  yields  the  angular 

^5  r\  Jl 

error . 

1-  =  (it)  (wSAT) 

=  (.7  sec)  (.05708  deg/sec)  =  .  03996  deg  (29) 

The  relationship  used  in  Equation  (24)  then  yields  the 
range  error. 

R  =  2  r  sin  (1  -  /  2  ) 

=  2  (7378  km)  sin  (.10998)  =  5.15  km  (30) 

While  this  particular  value  is  small,  it  does  partially 
account  for  the  difference  between  the  unperturbed  range 
and  the  perturbed  result. 

A  further  explanation  of  this  discrepancy  lies  in 
the  evolution  of  the  orbit's  eccentricity  over  time.  The 
time  variation  of  eccentricity  under  any  arbitrary  per¬ 
turbing  force  is  given  by  the  equation  (2:401) 
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1-e  sin  6 


a  (1-e  ) 


dt 


na 


2 

na  e 


-  r  i  r 


(31) 


/  3 

where  n  is  the  mean  motion,  vjj/a  ,  and  a  is  the  length 
of  the  semi-major  axis  of  the  orbit.  While  this  equation 
has  a  singularity  when  eccentricity  is  zero  and  therefore 
cannot  be  implemented  to  use  with  circular  orbits,  it  does 
show  that  if  there  are  any  force  components  in  the  local 
x-y  plane  then  the  eccentricity  of  the  orbit  will  change 
over  time.  For  any  non-circular  orbit  one  would  then 
expect  the  relative  position  to  change  with  time,  even 
for  unperturbed  motion.  Values  of  eccentricity  calculated 
during  the  baseline  case  run  were  on  the  order  of  10  3  to 
10'4. 


To  determine  how  much  of  the  range  oscillation 
shown  in  Figure  7  is  attributed  to  the  change  in  the  orbit's 
eccentricity,  a  case  was  run  using  a  typical  value  of 
eccentricity  from  the  baseline  run.  The  initial  condi¬ 
tions  for  this  non-circular  case  were  the  same  as  the 
baseline  (r  =  7378  km,  i  =  60,  ~2  ~  i  =  36°)  except  for 
the  eccentricity  which  was  set  at  .001.  The  results  are 
shown  in  Figure  8. 

Comparison  of  Figures  7  and  8  shows  that  the  out¬ 
comes  have  the  same  trend  for  both  the  initially  circular, 
perturbed  case  and  the  slightly  eccentric,  unperturbed 
case.  From  this  comparison  one  can  conclude  that  the 
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Fig.  8.  Unperturbed  Motion  with  Slightly  Eccentric  Orbits 

change  in  range  between  two  satellites  orbiting  within  the 
same  orbit  is  a  result  of  the  oscillation  in  the  orbit's 


eccentricity  caused  by  the  Earth's  oblateness  perturbation 
Examination  of  Figure  7  also  reveals  a  tendency 
for  the  range  oscillations  to  be  nonsymmetric  about  the 
average  value  rather  than  being  purely  sinusoidal. 

Equation  (31)  shows  that  the  time  rate  of  change  in  the 
orbit's  eccentricity  is  a  function  of  the  satellite’s  true 
anomally,  8.  However,  true  anomally  is  based  upon  the 
argument  of  perigee,  which  will  also  change  with  time  when 
the  orbit  is  perturbed  as  shown  in  the  next  equation 
(2:405)  . 


r  cot  i  sin  ( -+  ~) 


(32) 


F 

2  2  2 
na  v  l-e 

As  with  Equation  (31) ,  this  equation  is  unstable 
when  eccentricity  is  zero,  but  serves  to  show  that  the 
argument  of  perigee  is  influenced  by  local  force  components 
and  will  change  over  time.  The  result  of  gravitational 
perturbations  on  the  argument  of  perigee  is  a  rotation  of 
the  line  of  apsides.  The  magnitude  and  direction  of  this 
rotation  depends  upon  the  orbit's  radius  and  inclination 
(2:158).  Because  of  the  apsidal  rotation,  it  was  hypo¬ 
thesized  that  the  location  of  the  asymmetry  in  the  plot 
of  range  versus  time  was  dependent  upon  where  the  argument 
of  perigee  was  defined  in  the  initial  conditions. 

To  test  this  hypothesis  the  baseline  case  was  run 
several  times  using  different  initial  values  for  the 
argument  of  perigee,  but  with  the  satellite  always  at  the 
same  starting  point  in  the  orbit.  However,  the  results 
from  all  of  these  test  cases  yielded  the  same  plot  for 
range  versus  time,  leading  to  the  conclusion  that  the 
asymmetry  in  the  range  versus  time  plot  was  independent 
of  the  argument  of  perigee.  Therefore,  several  additional 
test  cases  were  formulated  that  used  a  common  argument  of 
perigee  but  started  with  the  satellites  at  various  locations 
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throughout  the  orbit.  Calculation  of  where  the  asymmetry 
and  the  maximum /minimum  range  values  occurred  relative  to 
the  ascending  node  showed  no  consistent  location  for  any  of 
these  three  conditions. 

The  only  definite  conclusions  which  can  be  drawn 
from  this  series  of  test  cases  are  that  the  range  between 
two  satellites  located  in  the  same  orbit  is  oscillatory, 
the  oscillations  repeat  with  a  frequency  equal  to  one 
orbital  period,  and  the  magnitude  of  the  oscillations 
neither  damps  out  nor  diverges  with  time. 

To  analyze  the  sensitivity  of  the  baseline  results 
to  variations  in  orbital  parameters,  several  test  cases 
were  run  using  different  combinations  of  orbital  radius 
and  inclination.  The  initial  conditions  for  all  of  the 
first  scenario  orbits  are  summarized  in  Table  I. 

The  results  for  all  of  these  cases  shewed  the  same 
general  trend  as  the  baseline  case  discussed  previously. 
Because  the  average  value  for  the  range  between  the  satel¬ 
lites  was  not  equal  to  the  calculated  range  for  a  separa¬ 
tion  of  36  degrees,  the  plots  of  range  versus  time  are 
difficult  to  compare  directly.  Therefore,  to  see  the  magni 
tude  of  the  oscillations  only,  the  results  were  normalized 
by  subtracting  each  test  case's  minimum  range  from  the 
calculated  range  at  each  time  step  during  that  run.  These 
normalized  results  are  shown  in  Figures  9  through  11, 
where  squares  indicate  the  40°  inclination  cases,  octagons 


TABLE  I 


INITIAL  ORBITAL 

ELEMENTS, 

FIRST 

SCENARIO 

Case 

£ 

e 

i 

ft 

CO 

1 

7378  km 

0 

O 

O 

rr 

45° 

0 

2* 

7378  km 

0 

60° 

45° 

0' 

3 

7378  km 

0 

00 

o 

o 

45° 

0' 

4 

6678  km 

0 

40° 

45° 

O' 

5 

6678  km 

0 

60° 

45° 

O' 

6 

6678  km 

0 

O 

O 

00 

45° 

o' 

7 

12378  km 

0 

o 

o 

O 

in 

o' 

8 

12378  km 

0 

60° 

45° 

O' 

0 


9 


12378  km 


0 


80° 


45 


indicate  60°  cases  and  triangles  indicate  80°  cases.  The 
total  magnitude  of  the  oscillations  ranged  from  2.5  km  to 
6.1  km  which,  for  all  cases,  is  a  deviation  of  less  than 
two-tenths  of  1  percent  of  the  average  value. 

Second  Scenario 

The  second  scenario  involved  evaluating  the  rela¬ 
tive  motion  between  two  satellites  in  similar,  but  differ¬ 
ent,  orbits  with  ascending  nodes  spaced  36  degrees  apart. 
Again,  this  spacing  equates  to  the  ten  orbit  constellation 
originally  proposed  as  the  baseline  condition.  Initial 
conditions  for  the  nine  test  cases  run  for  the  second 
scenario  are  outlined  in  Table  II,  with  Case  11  using  the 
baseline  orbit  conditions. 

TABLE  II 

INITIAL  ORBITAL  ELEMENTS,  SECOND  SCENARIO 


7378  km 
7378  km 
7378  km 
6678  km 
6678  km 
6678  km 
12378  km 
12378  km 
12378  km 
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Fig.  12.  Illustration  of  Second  Scenario 
Initial  Conditions 

Unlike  the  in-plane  arrangements  evaluated  under 
the  first  scenario,  the  range  between  satellites  in  differ¬ 
ent  orbit  planes  will  oscillate  between  some  maximum  and 
minimum  value  even  for  unperturbed  motion.  This  can  be 
illustrated  by  a  simple  example. 

Suppose  two  satellites  are  initially  positioned 
over  the  equator  (at  the  ascending  node)  in  circular  orbit 
that  are  both  at  90  degrees  of  inclination.  If  the  nodes 
of  the  two  orbit  planes  are  separated  by  some  amount,  the 
initial  range  can  be  visualized  as  illustrated  in  Figure  12 


The  range  at  this  point  can  be  calculated  as 

R  =  2  r  sin  (AQ/2) 

=  2  (7378)  sin ( 18 )  =  4559  km  (33) 

As  the  satellites  progress  in  their  orbits  from 
this  initial  position,  they  will  arrive  over  the  pole 
simultaneously  after  one-fourth  of  a  period,  at  which  time 
the  range  between  them  will  be  zero.  From  this  point  the 
range  will  again  increase  to  a  maximum  as  the  satellites 
reach  the  equator  at  the  descending  node.  This  pattern 
will  repeat  every  half  period.  For  orbits  at  inclinations 
other  than  90  degrees,  the  same  general  pattern  will  occur 
for  range  versus  time  but  the  minimum  range  value  will  be 
greater  than  zero  and  dependent  upon  the  inclination. 

For  the  baseline  orbits  inclined  at  60  degrees 
with  a  radius  of  7378  km,  the  minimum  range  was  calculated 
by  finding  the  inertial  position  vectors  for  the  two  satel¬ 
lites  from  the  orbital  elements  after  one  quarter  period. 
This  is  done  by  transforming  each  position  vector  expressed 
in  its  own  perifocal  reference  frame  into  the  inertial 
frame  (2:82,83)  and  then  finding  the  magnitude  of  the 
relative  position  vector.  This  calculation  resulted  in  an 
unperturbed  minimum  range  of  2279.972  km  for  the  baseline 


The  results  for  the  second  scenario  baseline  case 
are  shown  in  Figure  13.  The  maximum  range  for  this  case 
was  4559.543  km  and  the  minimum  was  2279.730  km. 

Figure  14  is  a  plot  of  range  versus  time  for  the 
unperturbed,  baseline  case.  The  maximum  range  for  this 
case  was  2279.973  km  and  the  minimum  was  4559.932  km. 
Comparison  with  Figure  13  shows  that  both  cases  are  stable 
over  time,  with  a  difference  between  maximum  and  minimum 
range  of  2279.959  km  for  the  unperturbed  case  and 
2279.813  km  with  perturbations. 

Because  of  the  large  excursions  in  range  present 
even  in  the  absence  of  perturbations,  the  results  are 
somewhat  difficult  to  interpret.  Because  the  range  rate 
for  the  second  scenario  cases  was  always  on  the  order  of 
hundreds  or  thousands  of  kilometers  per  second,  it  was 
not  possible  to  synchronize  the  results  of  the  perturbed 
and  unperturbed  cases  to  allow  the  deviations  to  be  calcu¬ 
lated.  However,  a  comparison  of  the  range  versus  time 
plots  for  the ■ perturbed  and  unperturbed  cases  shows  that 
the  difference  between  these  two  cases  is  insignificant. 

As  was  done  with  the  first  scenario  cases,  the 
baseline  case  for  the  second  scenario  was  tested  for  sensi¬ 
tivity  to  changes  in  orbital  altitude  and  inclination  using 
the  initial  conditions  shown  in  Table  II. 

The  results  of  these  additional  eight  cases  are 
shown  in  Figures  15  through  17,  and  as  with  the  baseline 
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Fig.  17.  Range  vs.  Time;  Cases  16,  17,  and  18 

case  the  range  oscillations  show  no  tendency  to  change 
with  time.  Again,  squares  correspond  to  40°  inclination 
orbits,  octagons  indicate  60°  orbits  and  triangles  indi¬ 
cate  80°  orbits.  This  indicates  that  the  relative  motion 


is  stable,  and,  as  calculated  for  the  baseline  condition, 
the  deviations  from  two-body  motion  are  minor. 


IV.  Conclusions  and  Recommendations 


Conclusion  s 

The  basic  question  this  research  effort  attempted 
to  answer  was  whether  the  perturbations  caused  by  the 
Earth's  non-uniform  gravitational  field  on  the  relative 
positions  of  satellites  arranged  as  a  constellation  are 
large  enough  to  significantly  alter  the  constellation's 
operational  capabilities.  Relative  motion  was  calculated 
using  a  set  of  equations  for  the  variation  of  the  differ¬ 
ence  between  the  orbital  elements  of  two  satellites.  These 
equations  allow  the  inclusion  of  any  arbitrary  perturbing 
forces,  and  such  forces  are  introduced  into  the  equations 
as  components  of  the  force  vector  expressed  in  a  reference 
frame  centered  at  the  satellite.  The  perturbing  force 
used  for  this  study  was  the  J2 harmonic  term  of  the  Earth's 
geopotential  function. 

Two  basic  scenarios  were  evaluated  for  the  study. 
The  first  scenario  looked  at  the  relative  motion  between 
two  satellites  orbiting  in  the  same  orbit  plane,  and  the 
second  scenario  calculated  relative  motion  between  satel¬ 
lites  in  similar  orbits  but  different  orbital  planes. 

The  results  presented  in  the  previous  chapter  show 
only  small  oscillations  in  range  between  adjacent  satel¬ 
lites  within  one  orbit  or  between  satellites  in  adjacent 


Because  of  similar  triangles,  can  be  found  from 


the  ratio 


R  2  r  sin  (A/2) 

P.  —  _ § _ 

R  2  r  sin  (A/2) 

s  s 


:  34 ) 


so  that 


R  =  R  — 
P  s  r 


35) 


Similarly,  oscillations  in  range  between  satellites  will 
result  in  oscillations  in  the  range  between  the  subpoints. 


,  R  =  A  R  — 
p  s  r 


36) 


Since  rg  will  always  be  greater  than  re,  ARp  will  always 

be  less  than  AR  .  When  combined  with  the  fact  that  the 
s 

range  oscillations  are  in  no  case  larger  than  4  km,  the 
effect  of  perturbation  induced  relative  motion  on  a  satel¬ 
lite  constellation  appears  to  be  insignificant.  The 
designer  of  a  satellite  constellation  should  be  able  to 
handle  motion  of  adjacent  weapons  or  sensor  footprints  by 
using  only  small  amounts  of  overlap  in  the  coverages. 

Keeping  in  mind  that  the  geopotential  model  used 
in  this  study  was  only  latitude  dependent,  the  findings  of 
this  investigation  point  to  the  conclusion  that  perturba¬ 
tions  caused  by  the  Earth's  non-uniform  gravitational  field 
on  the  relative  positions  of  satellites  arranged  as  a 


constellation  are  insignificant  and  do  not  need  to  be  con 
sidered  as  a  constraint  on  the  design  of  such  constella¬ 
tions. 

Recommendation s  for  Further  Study 

While  the  3^  harmonic  used  in  the  geopotential 
model  has  the  largest  magnitude  of  any  of  the  harmonic 
coefficients,  resonance  caused  by  longitudinal  variations 
in  the  earth's  gravity  field  are  possible,  but  no  longi¬ 
tudinal  dependence  was  modeled  in  this  study.  Therefore, 
expansion  of  the  geopotential  model  to  include  major 
sectoral  and/or  tesseral  harmonics  provides  an  avenue  for 
further  investigation  of  this  problem  which  may  have  sig- 
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Appendix:  Geopotential  Expansion  Theory 

The  dominant  cause  of  gravitational  perturbations 
on  earth  satellites  is  the  bulge  at  the  Earth's  equator. 
This  physical  anomally  shows  up  in  the  geopotential  expan¬ 
sion  of  Equation  (9)  as  the  term  for  n=2,  m=0  ,  and  can 
be  written  as 

V  =  GMR2P_°(cos  .  )  C _  A  / r  3  (A-l) 

e  2  20 

since  sin(m;)  vanishes  and  cos (m : )  equals  one  for  m=0 . 

m 

The  terms  P '  (cos  .  )  are  known  as  the  Legendre 
n 

polynomials,  and  for  this  case  can  be  solved  so  that 

P2°  (cos  )  =  (3  cos2  -1//2  (A- 2 ) 

Substituting  this  into  Equation  (.A-l)  results  m 

the  potential  being  a  function  of  latitude  only,  and  terms 

of  this  type  are  known  as  zonal  harmonics. 

The  harmonic  coefficients  C  and  S  are  derived 

r.m  nm 

empirically  from  satellite  observations  with  coefficients 
of  the  form  CnQ  usually  written  as  .  The  reported  value 
of  J2  is  1.0827x10  3  (9:290)  with  the  values  of  all  the 
other  harmonic  coefficients  being  on  the  order  of  10  ^ 
or  less.  Thus  the  J2  harmonic  has  the  largest  effect  on 
gravitational  perturbations. 
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